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Using highly efficient GPU-based simulations of the tight-binding Bogoliubov-de Gennes equa- 
tions we solve self-consistently for the pair correlation in rhombohedral (ABC) and Bernal (ABA) 
multilayer graphene by considering a finite intrinsic s-wave pairing potential. We find that the 
two different stacking configurations have opposite bulk/surface behavior for the order parameter. 
Surface superconductivity is robust for ABC stacked multilayer graphene even at very low pairing 
potentials for which the bulk order parameter vanishes, in agreement with a recent analytical ap- 
proach. In contrast, for Bernal stacked multilayer graphene, we find that the order parameter is 
always suppressed at the surface and that there exists a critical value for the pairing potential below 
which no superconducting order is achieved. We considered different doping scenarios and find that 
homogeneous doping strongly suppresses surface superconductivity while non-homogeneous field- 
induced doping has a much weaker effect on the superconducting order parameter. For multilayer 
structures with hybrid stacking (ABC and ABA) we find that when the thickness of each region is 
small (few layers), high-temperature surface superconductivity survives throughout the bulk due to 
the proximity effect between ABC/ABA interfaces where the order parameter is enhanced. 

PACS numbers: 73.43.-f, 73.23.-b. 73.63.-b 



I. INTRODUCTION 

Superconducting correlations in graphene-based struc- 
tures have been the focus of intensive theoretical and ex- 
perimental research even before graphene became one of 
the most important topics in condensed matter physics. 
Following last decade experimental evidences report- 
ing hints of superconductivity behavior in graphite^ 
and graphite intercalated compounda^r^, a considerable 
amount of theoretical studies have been devoted to pro- 
vide a clear understanding about possible mechanisms 
that could induce intrinsic superconducting states in sin- 
gle and multilayer graphene— More recently exper- 
imental investigations have reported intriguing traces 
of high-temperature superconducting behavior in highly 
oriented pyrolytic graphite (HOPG) samples^"—, feed- 
ing speculations about the existence of intrinsic super- 
conducting correlations in graphite and graphite-based 
compounds. Despite the fact that most of these ex- 
perimental evidence suggests that superconductivity in 
graphite compounds appears due to external causes, sev- 
eral theoretical studies reveal the possibility of inducing 
superconductivity in graphite by considering unconven- 
tional symmetry of the order parameter^^^'^"^. However, 
these calculations show that superconductivity becomes 
stable after considering disorder^ or high doping in the 
graphene layers^ while surface superconductivity ap- 
pears to be stable in clean rhombohedral graphite in the 
absence of external dopingii. Considering that most of 
these calculations are based on a reduced Hamiltonian or 
are performed within two-dimensional models, by ignor- 
ing the interplanar hopping, a numerical description of 
the superconducting correlation in multilayer graphene 
is urgently needed. 

In view of this, we provide in the following a numer- 



ical description of intrinsic superconductivity in mul- 
tilayer graphene at the tight-binding level. Following 
Ref. ml we consider a simple s-wave pairing symme- 
try in a ABC (or rhombohedral) stacking multilayered 
graphene structure. Calculations are also performed for 
ABA (or Bernal) stacking where its quadratic low-energy 
band structure^^ shows a remarkable difference from the 
\p\'^ momentum dependent band structure seen in ABC 
multilayer— structures. In particular we are interested in 
the limit of large number of layers (TV) where the lower- 
energy band in the rhombohedral case is flat over a large 
region in k-space signaling the suppression of the kinetic 
energy and therefore resulting in strong effect from in- 
teractions. In addition, because of the sensitive stack- 
ing dependent band structure in multilayer graphene, we 
also considered hybrid stacking cases. It is known that 
exfoliated few-layer graphite samples are usually found 
to exhibit very stable Bernal stacking but often also dis- 
play rhombohedral structures in part of the sampl o^^'^" . 
It is worth mentioning that during the preparation of the 
manuscript, new experimental results revealed the exis- 
tence of superconducting correlations at two-dimensional 
interfaces that appear when angle misalignments about 
the c-axis exist in HOPG-. 

By using highly efficient graphics card (GPU) simu- 
lations of the tight-binding Bogoliubov-de-Gennes equa- 
tions, we are able to solve self-consistently for the pair 
correlation in multilayer graphene by considering both 
planar and inter-planar coupling between nearest neigh- 
bors. Translational invariance is assumed along the 2- 
dimensional direction within the graphene sheets. In this 
way, an adequate description for the profile of supercon- 
ducting correlations along the direction perpendicular to 
the graphene layers is achieved. 

Our results confirm the main features of recent analyti- 



2 



cal approaches for ABC rhombohedral graphite where an 
enhancement of surface superconductivity, with respect 
to its bulk analog, was recently predicteclii. We find that 
the opposite behavior is present in ABA stacking where 
bulk dominates over surface superconductivity. This fact 
requires that a strong pairing potential is needed in order 
to obtain a non-zero pair correlation for the ABA case. 
In contrast, in the ABC lower pairing potential 

is sufficient to induce a large pair amplitude at the sur- 
face. In addition, we show how doping influences surface 
superconductivity, i.e. it strongly suppresses it. 

The paper is organized as follows. In Sec.|ll]we intro- 
duce the model and the numerical approach that we use. 
In Sec. Uni we present and discuss the results of our nu- 
merical calculations. Finally, we briefly summarize our 
flndings in Sec. |lVl 



II. MODEL AND CALCULATION APPROACH 

Superconducting correlations in multilayered graphene 
structures are described by the following mean-field sin- 
gle particle Hamiltonian written in Nambu space: 
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where the summation, is done over nearest neigh- 

bors within each layer while the summation, (/,m), is 
done for adjacent layers. The non-diagonal elements 
correspond to the s-wave superconducting order param- 
eter at the atomic site i in layer I while the diagonal 
elements "H;^ are the normal state components of the 
Hamiltonian: 
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FIG. 1. Schematic layout of the side view of a multilayer 
graphene structure stacked in two different configurations: 
Rhombohedral or ABC (left) and Bernal or ABA (right). 
Sublattices A and B are coupled by to within the same layer, 
while interlayer coupling is described by t. Integer coordinate 
z correspond to the index layer. 



where /i/ is the chemical potential and, according to the 
layout of Fig. [U nearest-neighbors sublattices A and B 
are coupled within the layers by the hopping parame- 
ter to «2.8eV while t = O.lio describes the hopping 
which couples A sites with the nearest B sites in the 
adjacent (upper and lower) layers. Bernal and rhom- 
bohedral stacking are defined according to the vertical 
symmetry along the z-axis as shown in Fig. [TJ Due 
to this symmetry, rhombohedral stacking allows us to 
reduce the description of the superconducting parame- 
ter to one of the sublattices, whereas the other sublat- 
tice can be deduced from a mirror refiection transforma- 
tion. For Bernal stacking, we consider a more practical 
way of sorting sites inside the supercell as dimmer sites, 
which correspond to the sublattices coupled by the in- 
terlayer hopping i, and no-dimmer sites or sublattices 
which does not participate in the coupling between ad- 
jacent graphene sheets. As we previously pointed out in 
a recent work^^, the pair correlation behaves differently 
in these inequivalent sites because at dimmer sites the 
density of states vanishes around the Dirac point while 
being finite at no-dimmer sites where the formation of 
Andreev states is more feasible. 

We will not specify the origin of superconductivity 
in the multilayer graphene structure, but rather assume 
= U{c\^c\t^ to be a conventional s-wave symmetric 
order paramelier and the pairing potential U is fixed and 
homogeneous in the whole structure. Under this assump- 
tion and considering translational invariance along the 
transversal directions, we solve self-consistently for the 
amplitude of the pair correlation | {c\^c\^ \ for the sub- 
lattice A (or both in the case of N-Bernal stacked lay- 
ers for N odd), in the z direction. We have considered 
graphene multi-layer supercells of size 42nmx25nmxA^ 
such that the order parameter is converged and no ad- 
ditional momentum summations in the parallel direc- 
tion are needed. Instead of a direct diagonalization of 
the Hamiltonian we performed the self-consistent mean- 
field calculation through a numerical approximation of 
the Gorkov Green's function by using the Chebyshev- 
Bogoliubov-de-Gennes metho d^^'^^ . Both, the normal 
and anomalous Gorkov Green's function, can be approxi- 
mated by a superposition of a finite number of Chebyshev 
polynomials as follows: 
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where the expansion coefficients for the diagonal, or nor- 
mal {a — 1), and the off-diagonal, or anomalous (a = 2), 
components of the 2x2 Green function are defined re- 
spectively as^^: 
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where Tn{x) = arccos[ncos(a;)] is the Chebyshev poly- 
nomial of order n, which satisfies the following recur- 
rence relation: T„+i(x) = 2xT{x) — r„_i(x). Once the 
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Hamiltonian has been normalized according toH H = 
{H — lb) /a, where the rescahng factors are a — {Emax — 
Eynin)/ (2 - ?y) and b = {Emax + £^m-i«)/2, with 77 > be- 
ing a smaU number, the expansion coefficients can be 
obtained through the dot product aljl^{n) - 

where (a| are the vectors (1| — (c;^| and (2| 
and the recursive vector = 2H\vn^i} — |vn-2)- We 
refer the reader to Ref. 2^ for more details about the 
numerical procedure. Since this iterative procedure in- 
volves mainly a successive application of the Hamiltonian 
matrix ([1]) on iterative vectors \vn), most of the compu- 
tational effort corresponds to sparse matrix- vector mul- 
tiplication which can be performed with high-efficiency 
by implementing parallel computations on GPUs by us- 
ing CUDA Nvidia. We are therefore able to solve effi- 
ciently multilayered graphene structures containing typ- 
ically several hundreds of thousand of atoms, which is 
half of the size of the Hamiltonian matrix ([1]). Physical 
quantities, like the local density of states and the pair 
correlation function, can be easily determined once the 



Green's functions are known: 
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III. NUMERICAL RESULTS 

We solved self-consistently for the order parameter 
along the z-direction, and show in Fig. [2][a) the different 
profiles of the order parameter (A^) for different values of 
the pairing potential {U < 0) in ABC stacked graphene 
with A^=20 layers. The order parameter is only shown for 
the A sublattice, while the B sublattice value is achieved 
by a mirror refiexion symmetry along z, as can be de- 
duced from Fig. [T] Analog results are shown in the inset 
of Fig. Ilja) for ABA stacking and different values of the 
pairing potential. 

We notice in Fig. [H^a), that the superconducting or- 
der parameter at the outermost layers is larger than the 
vanishing pair correlation in the bulk for all the U- values 
considered here. The same surface enhancement is ob- 
served when we decrease the pairing potential such that 
the penetration of the superconducting order parameter 
into the bulk becomes strongly suppressed. In the limit 
of very low pairing potential good agreement could be 
found with the analytical result previously reported in 
Ref. M 

On the other hand, an opposite surface-bulk supercon- 
ducting ratio is found for the Bernal stacking configura- 
tion (see inset in Fig[5l^a)). Self-consistent calculations 
performed for ABA show that the bulk value of the or- 
der parameter is dominant while surface superconduct- 
ing correlations are suppressed. Also, we can observed a 
sublattice polarization in the profile where pair corre- 
lation is found to be higher in non-dimmer sites as com- 
pared to dimmer sites with an energy difference which 
becomes smaller as the pairing potential is decreased. 
We return to this issue in a later discussion about the 



FIG. 2. (Color online) (a) Order parameter profile (A2) 
along the z-direction perpendicular to the graphene sheets 
(see Fig. [1} for various values of the s-wave attractive pairing 
potential U in ABC stacked multilayer graphene with A*' = 20 
layers. The inset shows the corresponding profile for the 
ABA case where dimmer and non-dimmer sites follow difi'er- 
ent curves for higher values of U. [/-dependence of the max- 
imum (b) and minimum (c) value of A^ for both ABC and 
ABA stacking configuration, d = 0.335nm is the inter-layer 
distance. 



density of states. 

A direct comparison between the surface and the bulk 
value of the pair correlation is given in Figs. [Hb) and 
djc) for both stacking configurations. Fig. [H^b) shows 
the maximum values of the superconducting correlation 
for both ABC and ABA cases for different values of U. 
According to the profiles observed in Fig. ^a), the max- 
imum value of the order parameter, maxjA^}, corre- 
sponds to the surface for ABC stacking while for the 
ABA stacking the maximum value corresponds to the 
bulk non-dimmer sites. In contrast, the log-linear repre- 
sentation presented in Fig. [Jfc) shows the U-dependence 
of the minimum value of the superconducting order pa- 
rameter, minjA^}, which corresponds to the bulk and 
surface locations for ABC and ABA, respectively. 

Two different regimes can be inferred from Figs. [2]^b) 
and [2jc) , depending whether the pairing potential U is 
larger or smaller than a critical value, |C/c|/io ^ 2.14. 
This value is also very close to the critical pairing for the 
ABA stacking, below which superconducting correlations 
vanish. As |C/| decreases, but is still larger than \Uc\, the 
order parameter decays exponentially for both ABC and 
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FIG. 3. (Color online) (a) Order parameter profile (Az) for a 
rhombohedral multilayer graphene consisting of various total 
number of layers A'^=10, 12, 16 and 20. (b) Order parameter 
profile for N=10 and 20 for different values of the interlayer 
hopping t/to=0.1, 0.092, 0.084, 0.076 and 0.068, in decreasing 
order as this is indicated by the arrows, (c) Dependence of the 
surface pair correlation, max{Az}, as a function of t. Here 
we used U = — 1.76to. 



ABA stacking, as the main contribution conies from the 
bulk. When \U\ < \Uc\, for ABC stacking, the bulk or- 
der parameter vanishes exponentially but the surface one 
is still finite. Even more interesting is the fact that the 
surface order parameter decays only linearly and is non- 
zero for all the values of |J7| that we considered in the 
simulation, down to |C^|/io = 1.76 giving a surface order 
parameter Amax/to ~ 0.002. The self-consistent calcula- 
tion becomes increasingly intensive as the order parame- 
ter decreases since more Chebyshev moments are needed 
to resolve the Green's function at higher and higher res- 
olutions. 

The bulk behavior resembles the superconducting crit- 
ical point reported for graphene, where in the undoped 
case it was found that the critical temperature vanishes 
below a finite value of the s-wave pairing interaction^"*. In 
contrast, the \U\ dependence of the surface order param- 
eter in the ABC stacking configuration suggests that the 
surface states, which form a flat band with suppressed 
kinetic energy, are strongly influenced by any exponen- 
tially small interaction. 

Since the pair correlation at the surface is always en- 
hanced for ABC stacking and survives even for lower val- 
ues of the pairing potential, we will further only also 



report numerical results for this stacking configuration. 

The dependence of the order parameter on the total 
number of graphene layers, N, and the interlayer cou- 
pling, t, is shown in Fig.[3]for U=-1.76tQ. An asymptotic 
enhancement of surface superconductivity is observed in 
Fig. EJa) as the Fermi surface size, defined by the flat 
band localized at the surfaces, increases with the num- 
ber of layers. On the other hand, the decrease of the 
interlayer coupling leads to the suppression of the order 
parameter as this is seen in Fig. ^h). The evolution of 
the surface pair correlation as a function of t is shown 
in Fig. ^c) and indicates an almost linear suppression 
as a consequence of the linear dependence of the Fermi 
surface size on t}^ 

In order to provide a better understanding of the pe- 
culiar behavior of the order parameter profile observed 
in Fig. [2] for the rhombohedral case, we next present the 
local density of states (LDOS). Fig. H shows the LDOS 
in different layers for both sublattices, A and B, within 
a small energy interval near the Fermi energy. The left 
panel of Fig. S] represents the LDOS at the surface where 
the superconducting order parameter is enhanced for sub- 
lattice A. There, the normal state LDOS shows a sharp 
peak due to the existence of flat bands with dispersion 
E ~ IpI^- The corresponding wave function of these 
states is localized at the surface and only on sublat- 
tice A. A zoomed-out view of the LDOS, over a wider 
range of energies, is shown in the inset of Fig. |4l This 
is very different for the B sublattice where the density 
of states vanishes around the Fermi energy and the su- 
perconducting coherence peaks are not visible. Despite 
this, a non-zero solution for the order parameter is ob- 
tained for atomic sites belonging to this sublattice as we 
can see in Figs. [Ifa) and [3] This non-zero solution ap- 
pears as a consequence of the proximity effect between 
the intra-layer neighbors A and B sites at the surface. 
We also observe less pronounced coherence peaks appear- 
ing in the LDOS of sublattice A in the layer adjacent to 
the surface (see central panel of Fig. U) while the LDOS 
vanishes for both sublattices in the bulk (see right panel 
of Fig. U]). According to this behavior of the LDOS for 
rhombohedral stacking, we expect that superconducting 
correlation will be more stable on the surface and on few 
adjacent layers rather than in the bulk where the LDOS 
vanishes around the Fermi energy. 
With respect to the LDOS in the Bernal case, it is well 
known that a sublattice polarization appears around the 
Fermi energy. A finite density of states is found for non- 
dimmer sites while the density of states vanishes at the 
dimmer sites^^. Such a polarization allows a finite order 
parameter to be induced only by large s-wave pairing po- 
tentials and therefore bulk superconductivity will not be 
stable for values of U for which surface superconductiv- 
ity in ABC is still finite (see Fig. [Sfb)). In this way, the 
suppression of surface superconductivity seen in Fig. [2jc) 
for ABA appears as a consequence of the lower density 
of states for surface non-dimmer sites when compared to 
its bulk value. 



5 



SURFACE SURFACE+1 BULK n/to[x1o1]: 2.0 3.1 

Vto[x10 3-5 2.7 0.7 




-0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 

Energy/to Energy/to Energy/tQ Energy/tg Energy/tg Energy/to 



FIG. 4. (Color online) Local density of states (LDOS) 
showing the formation of the s-wave superconducting gap for 
U — —2.0to at both sublattices A and B in different layers 
around the Fermi energy. Left and center panel show the 
LDOS at the surface and its adjacent layer while right panels 
show the LDOS at the bulk. The dashed line represent the 
corresponding normal state LDOS which shows the localized 
flat band at the outermost layers. The inset of the central 
panel shows the surface LDOS in ABC stacked graphene over 
a wider range of energies. 
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FIG. 5. (Color online) Effect of homogeneous doping on 
the intrinsic s-wave superconducting order parameter in a 
multilayer rhombohedral graphene with U = — 2.0to- We 
shown here half of the profile for different doping values, 
fi/to X 10^^=0, 0.8, 1.6, 2.4 and 3.1, in decreasing order indi- 
cates by the arrow. Inset: Evolution of surface superconduc- 
tivity (max{Az}) as a function of doping (/i) showing that 
an increase in fi leads to a strong suppression of the order 
parameter at the surface. 

In order to see the effect of external factors we have 
considered homogeneous doping, as well as an inho- 
mogeneous field-induced charge distribution, along the 
z-direction. The first case corresponds to graphite- 
intercalated compounds where dopant atoms are placed 
between the graphene layers while the inhomogeneous 



FIG. 6. (Color online) LDOS at the surface, on the A sub- 
lattice, for U = —2.0to and three different dopings considered 
before in Fig. [S] The normal state represented by the dashed 
line has been included for the corresponding values of doping 
fj,. As corresponds to maxjA^}. 

case can be easily realized in an experimental set-up 
where top and bottom electrodes have opposite gate 
voltage a^^i^^ . Fig. [5] shows the evolution of the order 
parameter profile as a function of the homogeneous dop- 
ing fiz = fJ- for a fixed value of the pairing potential 
U — — 2.0io- Notice that, surface superconductivity be- 
comes strongly suppressed as the doping shifts the Dirac 
point away from the Fermi energy. Total suppression oc- 
curs for doping lower than the value of the order param- 
eter at the surface. Looking at the LDOS we found that 
this critical doping coincides with the extinction of the 
coherence peaks and the rising of the single peak away 
form the Fermi energy. Despite the fact that we do not 
find any relevant difference between homogeneous doping 
and surface only doping in our self-consistent calculation, 
the value of the critical doping is slightly higher than the 
one reported by the analytical results where the critical 
doping was found to be ^crit=(2/3)A5. 

In addition to the homogeneous case we have also con- 
sidered inhomogeneous doping, achieved when an elec- 
tric field is applied perpendicular to the graphene sheets. 
Following Ref. [13 where the potential distribution due 
to an electric field was self-consistently calculated by 
taking into account screening effects, we consider the z- 
dependent doping as described in the inset of Fig. [71 For 
comparison purposes we have considered the same pair- 
ing potential U = — 2.0to, as we did in Fig. [Sj and three 
different cases for which the doping at the surfaces is 
strongly suppressing the superconducting correlations in 
the homogeneous case. Surprisingly, in contrast to homo- 
geneous doping, we found that the field-induced doping 
suppresses only slightly the pair correlation at the sur- 
face. According to this result, the electric-field induced 
gap in the inhomogeneous case is much lower than the 
value of the doping at the surface^^. Therefore, even 
considering the same surface doping in both cases, su- 
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FIG. 7. (Color online) Effect of inhomogeneous field-induced 
doping along the z-direction on the order parameter in a 
rhombohedral multilayer graphene. Some cases previously 
shown in Fig. [5] are included for comparison. Inset: Doping 
profile along the z-direction perpendicular to the graphene 
sheets. 
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FIG. 8. (Color online) LDOS showing the effect of field- 
induced doping along the z-direction (see inset in Fig. [7]) on 
surface superconductivity, obtained for U = — 2.0to in ABC 
multilayer graphene. Dashed curves represent the correspond- 
ing normal state LDOS. 



FIG. 9. (Color online) Self-consistent order parameter profile 
for a multilayer graphene considering the stacking configura- 
tion shown in the top part. Different point types correspond 
to different values of the pairing potential, U/to = 2.16, 2.12 
and 2.08, where U decrease in the direction pointed by the 
arrow. Open points represent the corresponding cases for an 
ordered ABC stacked multilayer graphene. 



ering stacking faults in multilayer graphene were recently 
reported^^, showing that surface superconductivity sur- 
vives at the interface between ABC and ABA stacking. 
However, superconductivity in the bulk is still seen be- 
ing suppressed when the thickness of the hybrid layers is 
large. By considering a few-layer structure where the ex- 
ternal layers have the ABC stacking configuration while 
the inner ones have ABA stacking, we expected a slight 
suppression of surface superconductivity and an enhance- 
ment of its bulk value. Fig. [9] shows self-consistent re- 
sults obtained for the hybrid structure depicted in the 
top part of the figure. While surface superconductivity 
appears suppressed as compared to the non-hybrid ABC 
multilayer, bulk correlation are also enhanced. On the 
basis of these results we suggest that this kind of hybrid 
stacked multilayer structure, or more complex combina- 
tions, could support high temperature superconductivity 
due to the interplay between surface superconductivity 
present in ABC stacking and the bulk superconductivity 
preserved in the ABA case. 



perconducting correlations are weakly affected by an in- 
homogeneous doping configuration. Fig. [8] shows how 
coherence peaks still survive for this doping configura- 
tion even if the value of the doping at the surface is on 
the order of the surface order parameter (as obtained for 
zero-doping). This is in contrast with the effect observed 
in Fig. [6] where for a similar level of doping causes a 
strong suppression of the superconducting gap. 

Finally, we consider multilayer graphene with hybrid 
stacking. In order to shift the high-temperature sur- 
face superconductivity, observed in ABC, to the bulk 
of the structure, we propose intercalated hybrid stacked 
graphite with few layers. Numerical calculations consid- 



IV. CONCLUSIONS 

By using a highly efficient CPU-based numerical pro- 
cedure we solved self-consistently for the s-wave order pa- 
rameter within a mean-field approach for a tight-binding 
model of ABA and ABC-stacked multilayer graphene. 
Main findings show that a surface superconducting state 
appear when the multilayer is in the ABC stacking config- 
uration. Opposite behavior is seen for the ABA stacking 
where bulk superconductivity is predominant but unsta- 
ble below a certain critical pairing potential. The LDOS 
for surface sites shows peculiar coherence peaks and sub- 
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lattice polarization, i.e. large LDOS in one sublattice and 
zero in the other. We showed that under homogeneous 
doping this state is quickly suppressed. We extracted 
a critical doping which is slightly higher than the one 
reported previously based on an approximate analytical 
studyii. In addition, we considered a field-induced inho- 
mogeneous doping and showed that the superconducting 
correlations survive in this case for higher values of iis- 
Finally we pointed out the importance of hybrid stack- 
ing structures where surface superconductivity related to 



ABC stacking could be preserved even in the bulk of the 
structure, suggesting a possible path for the survival of 
high temperature superconductivity. 
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